### TSCS Analysis ###

#Panelmatch version 2.0.1
#R version 4.2.0 
#Platform: x86_64-w64-mingw32/x64 (64-bit)

## Install/import Packages ##
library(foreign)
#install_github("insongkim/PanelMatch", dependencies=TRUE)
library(PanelMatch)



## NOTE: To replicate results with first episode only, replace MKactive with MKactive_alt

## Setting working directory, seed and clean R environment
rm(list=ls(all=TRUE))
set.seed(1234)

setwd("")

# Load data set and variable creation

# load data 
cohesion <- read.dta("datayearly.dta")
str(cohesion)

cohesion$dyadid <- as.integer(cohesion$dyadid)
cohesion$frag <- as.numeric(cohesion$frag)
cohesion$irregular <- as.numeric(cohesion$irregular2)
cohesion$rebsupnum <- as.numeric(cohesion$rebsupnum)
cohesion$excluded <- as.numeric(cohesion$excluded)
cohesion$prevactive<- as.numeric(cohesion$prevactive)
cohesion$fightcapnum<- as.numeric(cohesion$fightcapnum)
cohesion$terrcontnum<- as.numeric(cohesion$terrcontnum)
cohesion$conflictactive<- as.numeric(cohesion$conflictactive)
cohesion$year <- as.integer(cohesion$year)
cohesion$MKactive <- as.numeric(cohesion$MKactive)
cohesion$exclpopn <- as.numeric(cohesion$exclpopn)
cohesion$e_v2x_neopat <- as.numeric(cohesion$e_v2x_neopat)
cohesion$dyadsfirstyearconflict <- as.numeric(cohesion$dyadsfirstyearconflict)
cohesion$MKprevious <- as.numeric(cohesion$MKprevious)

class(cohesion$year)
colnames(cohesion) 

### Appendix Table 23 ###

#placeholder matrix for the results
attbyperiod <- matrix(NA, nrow=24, ncol=4)
colnames(attbyperiod) <- c("Estimate", "Std. Error", "2.5%", "97.5%")
rownames(attbyperiod) <- rep(c("t+1", "t+2", "t+3"),8)

## ATT by period 1 lag with missing values
PM.results1lagmain <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid", 
                                 treatment = "MKactive", refinement.method = "ps.weight", 
                                 data = cohesion, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(irregular, 1:1)) + I(lag(rebsupnum, 1:1)) + I(lag(excluded, 1:1)) + I(lag(prevactive, 1:1)) + I(lag(fightcapnum, 1:1)) + I(lag(terrcontnum, 1:1)) + I(lag(e_v2x_neopat, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(conflictactive, 1:1)),
                                 qoi = "att"
                                 , outcome.var = "frag",
                                 lead = 1:3, forbid.treatment.reversal = TRUE)


mset.obj1lagmain <- PM.results1lagmain$att
attributes(mset.obj1lagmain)
summary(mset.obj1lagmain)

# covariate balance
get_covariate_balance(PM.results1lagmain$att, cohesion, covariates = c("irregular", "rebsupnum", "excluded", "prevactive", "fightcapnum", "terrcontnum", "e_v2x_neopat", "exclpopn", "dyadsfirstyearconflict", "conflictactive", "frag"), verbose = T, plot = F,
                      reference.line = F, legend = TRUE, ylab = "Time")


# estimate effect
set.seed(1234)
PE.results1lagmain <- PanelEstimate(number.iterations = 1000, sets = PM.results1lagmain, 
                                    data = cohesion)

summary(PE.results1lagmain)

# replacing estimates in the placeholder matrix
attbyperiod[1:3,1] <- round(PE.results1lagmain$estimates, 3)
attbyperiod[1:3,2] <- round(PE.results1lagmain$standard.error, 3)
attbyperiod[1,3:4] <- round(quantile(PE.results1lagmain$bootstrapped.estimates[,1], c(0.025, 0.975)),3)
attbyperiod[2,3:4] <- round(quantile(PE.results1lagmain$bootstrapped.estimates[,2], c(0.025, 0.975)),3)
attbyperiod[3,3:4] <- round(quantile(PE.results1lagmain$bootstrapped.estimates[,3], c(0.025, 0.975)),3)


## ATT by period 1 lag with missing values
PM.results1lagmissingmain <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid", 
                                        treatment = "MKactive", refinement.method = "ps.weight", 
                                        data = cohesion, match.missing = TRUE, 
                                        covs.formula = ~ I(lag(irregular, 1:1)) + I(lag(rebsupnum, 1:1)) + I(lag(excluded, 1:1)) + I(lag(prevactive, 1:1)) + I(lag(fightcapnum, 1:1)) + I(lag(terrcontnum, 1:1)) + I(lag(e_v2x_neopat, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(conflictactive, 1:1)),
                                        qoi = "att"
                                        , outcome.var = "frag",
                                        lead = 1:3, forbid.treatment.reversal = TRUE)

mset.obj1lagmissmain <- PM.results1lagmissingmain$att
attributes(mset.obj1lagmissmain)
summary(mset.obj1lagmissmain)

#covariate balance

get_covariate_balance(PM.results1lagmissingmain$att, cohesion, covariates = c("irregular", "rebsupnum", "excluded", "prevactive", "fightcapnum", "terrcontnum", "e_v2x_neopat", "exclpopn", "dyadsfirstyearconflict", "conflictactive", "frag"), verbose = T, plot = F,
                      reference.line = F, legend = TRUE, ylab = "Standardized Mean Difference after Refinement", xlab = "Time")

#estimate effect
set.seed(1234)
PE.results1lagmissingmain <- PanelEstimate(number.iterations = 1000, sets = PM.results1lagmissingmain, 
                                           data = cohesion)
summary(PE.results1lagmissingmain)

# replacing estimates in the placeholder matrix
attbyperiod[4:6,1] <- round(PE.results1lagmissingmain$estimates, 3)
attbyperiod[4:6,2] <- round(PE.results1lagmissingmain$standard.error, 3)
attbyperiod[4,3:4] <- round(quantile(PE.results1lagmissingmain$bootstrapped.estimates[,1], c(0.025, 0.975)),3)
attbyperiod[5,3:4] <- round(quantile(PE.results1lagmissingmain$bootstrapped.estimates[,2], c(0.025, 0.975)),3)
attbyperiod[6,3:4] <- round(quantile(PE.results1lagmissingmain$bootstrapped.estimates[,3], c(0.025, 0.975)),3)


## ATT by period 1 lag, previous MKs

PM.results1lagmain <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid", 
                                 treatment = "MKactive", refinement.method = "ps.weight", 
                                 data = cohesion, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(irregular, 1:1)) + I(lag(rebsupnum, 1:1)) + I(lag(excluded, 1:1)) + I(lag(prevactive, 1:1)) + I(lag(fightcapnum, 1:1)) + I(lag(terrcontnum, 1:1)) + I(lag(e_v2x_neopat, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(conflictactive, 1:1)) + I(lag(MKprevious, 1:1)),
                                 qoi = "att"
                                 , outcome.var = "frag",
                                 lead = 1:3, forbid.treatment.reversal = TRUE)

mset.obj1lagmain <- PM.results1lagmain$att
attributes(mset.obj1lagmain)
summary(mset.obj1lagmain)

#covariate balance

get_covariate_balance(PM.results1lagmain$att, cohesion, covariates = c("irregular", "rebsupnum", "excluded", "prevactive", "fightcapnum", "terrcontnum", "e_v2x_neopat", "exclpopn", "dyadsfirstyearconflict", "conflictactive", "frag", "MKprevious"), verbose = T, plot = F,
                      reference.line = F, legend = TRUE, ylab = "Time")

#estimate effect
set.seed(1234)
PE.results1lagmain <- PanelEstimate(number.iterations = 1000, sets = PM.results1lagmain, 
                                    data = cohesion)
summary(PE.results1lagmain)

# replacing estimates in the placeholder matrix
attbyperiod[7:9,1] <- round(PE.results1lagmain$estimates, 3)
attbyperiod[7:9,2] <- round(PE.results1lagmain$standard.error, 3)
attbyperiod[7,3:4] <- round(quantile(PE.results1lagmain$bootstrapped.estimates[,1], c(0.025, 0.975)),3)
attbyperiod[8,3:4] <- round(quantile(PE.results1lagmain$bootstrapped.estimates[,2], c(0.025, 0.975)),3)
attbyperiod[9,3:4] <- round(quantile(PE.results1lagmain$bootstrapped.estimates[,3], c(0.025, 0.975)),3)

## ATT by period 1 lag, previous MKs with missing values

PM.results1lagmissingmain <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid", 
                                        treatment = "MKactive", refinement.method = "ps.weight", 
                                        data = cohesion, match.missing = TRUE, 
                                        covs.formula = ~ I(lag(irregular, 1:1)) + I(lag(rebsupnum, 1:1)) + I(lag(excluded, 1:1)) + I(lag(prevactive, 1:1)) + I(lag(fightcapnum, 1:1)) + I(lag(terrcontnum, 1:1)) + I(lag(e_v2x_neopat, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(conflictactive, 1:1)) + I(lag(MKprevious, 1:1)),
                                        qoi = "att"
                                        , outcome.var = "frag",
                                        lead = 1:3, forbid.treatment.reversal = TRUE)

mset.obj1lagmissmain <- PM.results1lagmissingmain$att
attributes(mset.obj1lagmissmain)
summary(mset.obj1lagmissmain)

#covariate balance

get_covariate_balance(PM.results1lagmissingmain$att, cohesion, covariates = c("irregular", "rebsupnum", "excluded", "prevactive", "fightcapnum", "terrcontnum", "e_v2x_neopat", "exclpopn", "dyadsfirstyearconflict", "conflictactive", "frag", "MKprevious"), verbose = T, plot = F,
                      reference.line = F, legend = TRUE, ylab = "Standardized Mean Difference after Refinement", xlab = "Time")


#estimate effect
set.seed(1234)
PE.results1lagmissingmain <- PanelEstimate(number.iterations = 1000, sets = PM.results1lagmissingmain, 
                                           data = cohesion)
summary(PE.results1lagmissingmain)

# replacing estimates in the placeholder matrix
attbyperiod[10:12,1] <- round(PE.results1lagmissingmain$estimates, 3)
attbyperiod[10:12,2] <- round(PE.results1lagmissingmain$standard.error, 3)
attbyperiod[10,3:4] <- round(quantile(PE.results1lagmissingmain$bootstrapped.estimates[,1], c(0.025, 0.975)),3)
attbyperiod[11,3:4] <- round(quantile(PE.results1lagmissingmain$bootstrapped.estimates[,2], c(0.025, 0.975)),3)
attbyperiod[12,3:4] <- round(quantile(PE.results1lagmissingmain$bootstrapped.estimates[,3], c(0.025, 0.975)),3)

## ATT by period 1 lag, past outcome

PM.results1lagmain <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid", 
                                 treatment = "MKactive", refinement.method = "ps.weight", 
                                 data = cohesion, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(irregular, 1:1)) + I(lag(rebsupnum, 1:1)) + I(lag(excluded, 1:1)) + I(lag(prevactive, 1:1)) + I(lag(fightcapnum, 1:1)) + I(lag(terrcontnum, 1:1)) + I(lag(e_v2x_neopat, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(conflictactive, 1:1)) + I(lag(frag, 1:1)),
                                 qoi = "att"
                                 , outcome.var = "frag", 
                                 lead = 1:3,  forbid.treatment.reversal = TRUE)

mset.obj1lagmain <- PM.results1lagmain$att
attributes(mset.obj1lagmain)
summary(mset.obj1lagmain)

#covariate balance

get_covariate_balance(PM.results1lagmain$att, cohesion, covariates = c("irregular", "rebsupnum", "excluded", "prevactive", "fightcapnum", "terrcontnum", "e_v2x_neopat", "exclpopn", "dyadsfirstyearconflict", "conflictactive", "frag"), verbose = T, plot = F,
                      reference.line = F, legend = TRUE, ylab = "Time")

#estimate effect
set.seed(1234)
PE.results1lagmain <- PanelEstimate(number.iterations = 1000, sets = PM.results1lagmain, 
                                    data = cohesion)
summary(PE.results1lagmain)

# replacing estimates in the placeholder matrix
attbyperiod[13:15,1] <- round(PE.results1lagmain$estimates, 3)
attbyperiod[13:15,2] <- round(PE.results1lagmain$standard.error, 3)
attbyperiod[13,3:4] <- round(quantile(PE.results1lagmain$bootstrapped.estimates[,1], c(0.025, 0.975)),3)
attbyperiod[14,3:4] <- round(quantile(PE.results1lagmain$bootstrapped.estimates[,2], c(0.025, 0.975)),3)
attbyperiod[15,3:4] <- round(quantile(PE.results1lagmain$bootstrapped.estimates[,3], c(0.025, 0.975)),3)

## ATT by period 1 lag, past outcome with missing values

PM.results1lagmissingmain <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid", 
                                        treatment = "MKactive", refinement.method = "ps.weight", 
                                        data = cohesion, match.missing = TRUE, 
                                        covs.formula = ~ I(lag(irregular, 1:1)) + I(lag(rebsupnum, 1:1)) + I(lag(excluded, 1:1)) + I(lag(prevactive, 1:1)) + I(lag(fightcapnum, 1:1)) + I(lag(terrcontnum, 1:1)) + I(lag(e_v2x_neopat, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(conflictactive, 1:1)) + I(lag(frag, 1:1)),
                                        qoi = "att"
                                        , outcome.var = "frag", 
                                        lead = 1:3,  forbid.treatment.reversal = TRUE)

mset.obj1lagmissmain <- PM.results1lagmissingmain$att
attributes(mset.obj1lagmissmain)
summary(mset.obj1lagmissmain)

#covariate balance

get_covariate_balance(PM.results1lagmissingmain$att, cohesion, covariates = c("irregular", "rebsupnum", "excluded", "prevactive", "fightcapnum", "terrcontnum", "e_v2x_neopat", "exclpopn", "dyadsfirstyearconflict", "conflictactive", "frag"), verbose = T, plot = F,
                      reference.line = F, legend = TRUE, ylab = "Standardized Mean Difference after Refinement", xlab = "Time")


#estimate effect
set.seed(1234)
PE.results1lagmissingmain <- PanelEstimate(number.iterations = 1000, sets = PM.results1lagmissingmain, 
                                           data = cohesion)
summary(PE.results1lagmissingmain)

# replacing estimates in the placeholder matrix
attbyperiod[16:18,1] <- round(PE.results1lagmissingmain$estimates, 3)
attbyperiod[16:18,2] <- round(PE.results1lagmissingmain$standard.error, 3)
attbyperiod[16,3:4] <- round(quantile(PE.results1lagmissingmain$bootstrapped.estimates[,1], c(0.025, 0.975)),3)
attbyperiod[17,3:4] <- round(quantile(PE.results1lagmissingmain$bootstrapped.estimates[,2], c(0.025, 0.975)),3)
attbyperiod[18,3:4] <- round(quantile(PE.results1lagmissingmain$bootstrapped.estimates[,3], c(0.025, 0.975)),3)

## ATT by period 1 lag, previous MKs, past outcome

PM.results1lagmain <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid", 
                                 treatment = "MKactive", refinement.method = "ps.weight", 
                                 data = cohesion, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(irregular, 1:1)) + I(lag(rebsupnum, 1:1)) + I(lag(excluded, 1:1)) + I(lag(prevactive, 1:1)) + I(lag(fightcapnum, 1:1)) + I(lag(terrcontnum, 1:1)) + I(lag(e_v2x_neopat, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(conflictactive, 1:1))  + I(lag(frag, 1:1)) + I(lag(MKprevious, 1:1)),
                                 qoi = "att"
                                 , outcome.var = "frag",
                                 lead = 1:3, forbid.treatment.reversal = TRUE)

mset.obj1lagmain <- PM.results1lagmain$att
attributes(mset.obj1lagmain)
summary(mset.obj1lagmain)

#Balance

get_covariate_balance(PM.results1lagmain$att, cohesion, covariates = c("irregular", "rebsupnum", "excluded", "prevactive", "fightcapnum", "terrcontnum", "e_v2x_neopat", "exclpopn", "dyadsfirstyearconflict", "conflictactive", "frag", "MKprevious"), verbose = T, plot = F,
                      reference.line = F, legend = TRUE, ylab = "Time")

#estimate effect
set.seed(1234)
PE.results1lagmain <- PanelEstimate(number.iterations = 1000, sets = PM.results1lagmain, 
                                    data = cohesion)
summary(PE.results1lagmain)


# replacing estimates in the placeholder matrix
attbyperiod[19:21,1] <- round(PE.results1lagmain$estimates, 3)
attbyperiod[19:21,2] <- round(PE.results1lagmain$standard.error, 3)
attbyperiod[19,3:4] <- round(quantile(PE.results1lagmain$bootstrapped.estimates[,1], c(0.025, 0.975)),3)
attbyperiod[20,3:4] <- round(quantile(PE.results1lagmain$bootstrapped.estimates[,2], c(0.025, 0.975)),3)
attbyperiod[21,3:4] <- round(quantile(PE.results1lagmain$bootstrapped.estimates[,3], c(0.025, 0.975)),3)

## ATT by period 1 lag, previous MKs, past outcome with missing values AND Figure 2

PM.results1lagmissingmain <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid", 
                                        treatment = "MKactive", refinement.method = "ps.weight", 
                                        data = cohesion, match.missing = TRUE, 
                                        covs.formula = ~ I(lag(irregular, 1:1)) + I(lag(rebsupnum, 1:1)) + I(lag(excluded, 1:1)) + I(lag(prevactive, 1:1)) + I(lag(fightcapnum, 1:1)) + I(lag(terrcontnum, 1:1)) + I(lag(e_v2x_neopat, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(conflictactive, 1:1))  + I(lag(frag, 1:1)) + I(lag(MKprevious, 1:1)),
                                        qoi = "att"
                                        , outcome.var = "frag",
                                        lead = 1:3, forbid.treatment.reversal = TRUE)

PM.results1lagmissingnone <- PanelMatch(lag = 1, time.id = "year", unit.id = "dyadid", 
                                        treatment = "MKactive", refinement.method = "none", 
                                        data = cohesion, match.missing = TRUE, 
                                        covs.formula = ~ I(lag(irregular, 1:1)) + I(lag(rebsupnum, 1:1)) + I(lag(excluded, 1:1)) + I(lag(prevactive, 1:1)) + I(lag(fightcapnum, 1:1)) + I(lag(terrcontnum, 1:1)) + I(lag(e_v2x_neopat, 1:1)) + I(lag(dyadsfirstyearconflict, 1:1)) + I(lag(exclpopn, 1:1)) + I(lag(conflictactive, 1:1))  + I(lag(frag, 1:1)) + I(lag(MKprevious, 1:1)),
                                        qoi = "att"
                                        , outcome.var = "frag",
                                        lead = 1:3, forbid.treatment.reversal = TRUE)


mset.obj1lagmissmain <- PM.results1lagmissingmain$att
attributes(mset.obj1lagmissmain)
summary(mset.obj1lagmissmain)

#covariate balance

get_covariate_balance(PM.results1lagmissingnone$att, cohesion, covariates = c("irregular", "rebsupnum", "excluded", "prevactive", "fightcapnum", "terrcontnum", "e_v2x_neopat", "exclpopn", "dyadsfirstyearconflict", "conflictactive", "frag", "MKprevious"), verbose = T, plot = F,
                      reference.line = F, legend = TRUE, ylab = "Standardized Mean Difference after Refinement", xlab = "Time")


get_covariate_balance(PM.results1lagmissingmain$att, cohesion, covariates = c("irregular", "rebsupnum", "excluded", "prevactive", "fightcapnum", "terrcontnum", "e_v2x_neopat", "exclpopn", "dyadsfirstyearconflict", "conflictactive", "frag", "MKprevious"), verbose = T, plot = F,
                      reference.line = F, legend = TRUE, ylab = "Standardized Mean Difference after Refinement", xlab = "Time")

#compare refined and non-refined set
balance_scatter(non_refined_set = PM.results1lagmissingnone$att,
                matched_set_list = list(PM.results1lagmissingmain$att),
                data = cohesion,
                covariates = c("irregular", "rebsupnum", "excluded","prevactive", "fightcapnum", "terrcontnum", "e_v2x_neopat", "dyadsfirstyearconflict", "exclpopn", "conflictactive", "frag", "MKprevious"))

### Figure 2

pdf("figure2.pdf")
scatter4b<-balance_scatter(non_refined_set = PM.results1lagmissingnone$att,
                           matched_set_list = list(PM.results1lagmissingmain$att),
                           data = cohesion,
                           covariates = c("irregular", "rebsupnum", "excluded","prevactive", "fightcapnum", "terrcontnum", "e_v2x_neopat", "dyadsfirstyearconflict", "exclpopn", "conflictactive", "frag", "MKprevious"))
dev.off()

# estimate effect
set.seed(1234)
PE.results1lagmissingmain <- PanelEstimate(number.iterations = 1000, sets = PM.results1lagmissingmain, 
                                           data = cohesion)
summary(PE.results1lagmissingmain)

# replacing estimates in the placeholder matrix
attbyperiod[22:24,1] <- round(PE.results1lagmissingmain$estimates, 3)
attbyperiod[22:24,2] <- round(PE.results1lagmissingmain$standard.error, 3)
attbyperiod[22,3:4] <- round(quantile(PE.results1lagmissingmain$bootstrapped.estimates[,1], c(0.025, 0.975)),3)
attbyperiod[23,3:4] <- round(quantile(PE.results1lagmissingmain$bootstrapped.estimates[,2], c(0.025, 0.975)),3)
attbyperiod[24,3:4] <- round(quantile(PE.results1lagmissingmain$bootstrapped.estimates[,3], c(0.025, 0.975)),3)


View(attbyperiod)

write.csv(attbyperiod, "ATT.csv")

